term <- volume/max_depth/pi*3-Radius^2
radius_imaginary <- polyroot(c(-term,Radius,1))
real<-Re(radius_imaginary)
radius <-real[which(real>0)]
if (length(radius)==0){
a=2
print("Error 2")
return(a)
}
stopifnot(length(radius)>0)
if (Radius<radius){
a=1
print("Error 1")
return(a)
}
stopifnot(Radius>=radius)
area_layer <- pi*Radius^2
hipsographic[1,1]=0
hipsographic[1,2]=area_layer
hipsographic[1,3]=volume
slope <- (Radius-radius)/max_depth
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 1:floor(max_depth)){
new_radius <- Radius-i*slope
area_layer <- pi*new_radius^2
hipsographic[i+1,1]=i
hipsographic[i+1,2]=area_layer
hipsographic[i+1,3]=pi/3*(max_depth-i)*(new_radius^2+new_radius*radius+radius^2)
}
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= pi*radius^2
hipsographic[i+2,3]= 0
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= pi*radius^2
hipsographic[2,3]= 0
}
}
a=0
return(a)
#return(hipsographic)
}
library(foreign)
data<-read.dbf("/home/rmarce/Cloud/a. WATExR/ISIMIP/Area_depth February 2022/Hydrolakes/HydroLAKES_polys_v10_shp/HydroLAKES_polys_v10.dbf")
# Function for the hypsographic curve assuming a truncated cone from volume (m^3), surface area (m^2), and maximum depth (m).
# Calculated every meter from surface and at maximum depth.
hipsographicR_TruncatedCone <-  function(volume,surface_area,max_depth){
#volume=total volume of lake (m^3)
#surface_area= surface area (m^2)
#max_depth=maximum depth (m)
#initializing variable
##fields for depth,area,accumulated volume
#conditional for including the maximum depth in case it is not an integer depth
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
#solving a truncated cone
Radius <- sqrt(surface_area/pi) #radius of the large end of the truncated cone (surface of the lake)
term <- volume/max_depth/pi*3-Radius^2 #solving the truncated cone equation [=0])
radius_imaginary <- polyroot(c(-term,Radius,1)) #finding the roots of the equation
real<-Re(radius_imaginary) # taking the real part of the result
#the positive real result is the result we want:
radius <-real[which(real>0)] #radius of the small end of the truncated cone (bottom of the lake)
#checks if the solution makes physical sense and print an error otherwise:
#the root must have a real part
if (length(radius)==0){
a=2
print("Error 2")
return(a)
}
stopifnot(length(radius)>0)
#surface radius must be larger or equal than bottom radius
if (Radius<radius){
a=1
print("Error 1")
return(a)
}
stopifnot(Radius>=radius)
#Surface area
area_layer <- pi*Radius^2 #must be equal to surface_area
#Result at the surface
hipsographic[1,1]=0
hipsographic[1,2]=area_layer
hipsographic[1,3]=volume
slope <- (Radius-radius)/max_depth #for calculating the radius at each depth
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 1:floor(max_depth)){
new_radius <- Radius-i*slope
area_layer <- pi*new_radius^2
hipsographic[i+1,1]=i
hipsographic[i+1,2]=area_layer
hipsographic[i+1,3]=pi/3*(max_depth-i)*(new_radius^2+new_radius*radius+radius^2)
}
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= pi*radius^2
hipsographic[i+2,3]= 0
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= pi*radius^2
hipsographic[2,3]= 0
}
}
#a=0
#return(a)
return(hipsographic)
}
volume <- 168555237 #m3
surface_area <- 5753968 #m2
max_depth <- 60.54 #m
# Function for the hypsographic curve assuming a truncated cone from volume (m^3), surface area (m^2), and maximum depth (m).
# Calculated every meter from surface and at maximum depth.
hipsographicR_TruncatedCone <-  function(volume,surface_area,max_depth){
#volume=total volume of lake (m^3)
#surface_area= surface area (m^2)
#max_depth=maximum depth (m)
#initializing variable
##fields for depth,area,accumulated volume
#conditional for including the maximum depth in case it is not an integer depth
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
#solving a truncated cone
Radius <- sqrt(surface_area/pi) #radius of the large end of the truncated cone (surface of the lake)
term <- volume/max_depth/pi*3-Radius^2 #solving the truncated cone equation [=0])
radius_imaginary <- polyroot(c(-term,Radius,1)) #finding the roots of the equation
real<-Re(radius_imaginary) # taking the real part of the result
#the positive real result is the result we want:
radius <-real[which(real>0)] #radius of the small end of the truncated cone (bottom of the lake)
#checks if the solution makes physical sense and print an error otherwise:
#the root must have a real part
if (length(radius)==0){
a=2
print("Error 2")
return(a)
}
stopifnot(length(radius)>0)
#surface radius must be larger or equal than bottom radius
if (Radius<radius){
a=1
print("Error 1")
return(a)
}
stopifnot(Radius>=radius)
#Surface area
area_layer <- pi*Radius^2 #must be equal to surface_area
#Result at the surface
hipsographic[1,1]=0
hipsographic[1,2]=area_layer
hipsographic[1,3]=volume
slope <- (Radius-radius)/max_depth #for calculating the radius at each depth
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 1:floor(max_depth)){
new_radius <- Radius-i*slope
area_layer <- pi*new_radius^2
hipsographic[i+1,1]=i
hipsographic[i+1,2]=area_layer
hipsographic[i+1,3]=pi/3*(max_depth-i)*(new_radius^2+new_radius*radius+radius^2)
}
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= pi*radius^2
hipsographic[i+2,3]= 0
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= pi*radius^2
hipsographic[2,3]= 0
}
}
#a=0
#return(a)
return(hipsographic)
}
hipsographicR_TruncatedCone(volume,surface_area,max_depth)
hipso_TC <-  hipsographicR_TruncatedCone(volume,surface_area,max_depth)
plot(hipso_TC[,2],hipso_TC[,1])
plot(hipso_TC[,2],hipso_TC[,1],type="l")
plot(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(y)))
plot(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,2])))
plot(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,1])))
z_m=max_depth*0.47
Vd <- 3*z_m/max_depth
b <- log(Vd)
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
mean_depth <- volume/surface
mean_depth <- volume/surface_area
z_m <- max_depth*0.47
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log(Vd)
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
i in 0:floor(max_depth)
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log10(Vd)
#if per les diferents possibilitats
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
Vd=1.86
b <- log10(Vd)
Hd <- 10^(-19261*b^5+27179*b^4-15506*b^3+4432:6*b^2-639.51*b+36.442)
Hd <- 10^(-19261*b^5+27179*b^4-15506*b^3+4432.6*b^2-639.51*b+36.442)
b <- log(Vd)
Hd <- 10^(-19261*b^5+27179*b^4-15506*b^3+4432.6*b^2-639.51*b+36.442)
b <- log10(Vd)
Hd <- 10^(-19261*b^5+27179*b^4-15506*b^3+4432.6*b^2-639.51*b+36.442)
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log10(Vd)
#if per les diferents possibilitats
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
i=0
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
View(hipsographic)
surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-Hd^-1)^2)) * (  (2*max_depth*Hd^(-i/max_depth-1))/log(Hd) - (max_depth*Hd^(-2*i/max_depth))/(2*log(Hd)) + (Hd^-2)*i )
View(hipsographic)
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
View(hipsographic)
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
i
(surface_area/((1-(Hd^-1))^2))
(surface_area/((1-(Hd^-1))^2))
(surface_area/((1-(Hd^-1))^2))
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
volume <- 168555237 #m3
surface_area <- 5753968 #m2
max_depth <- 60.54 #m
mean_depth <- volume/surface_area
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log10(Vd)
#if per les diferents possibilitats
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
i=0
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
}
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
}
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
}
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log10(Vd)
#if per les diferents possibilitats
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
}
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
}
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
lines(hipsographic[,2],hipsographic[,1], col="red")
plot(hipsographic[,2],hipsographic[,1], col="red",ylim = rev(range(hipsographic[,1])))
lines(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,1])))
sau <- c(0	5999766.8
0.54	5944169.15
1.54	5810168.49
2.54	5661639.35
3.54	5534433.79
4.54	5399721.46
5.54	5264874.32
6.54	5097163.86
7.54	4960082.47
8.54	4843177.19
9.54	4720941.89
10.54	4608242.26
11.54	4495187.73
12.54	4393014.88
13.54	4292236.4
14.54	4193774.53
15.54	4084796.54
16.54	3930219.88
17.54	3790900.21
18.54	3671706.4
19.54	3563206.57
20.54	3426626.36
21.54	3326417.75
22.54	3229475.01
23.54	3135233.67
24.54	3042020.04
25.54	2945708.42
26.54	2846967.51
27.54	2724520.41
28.54	2603570.27
29.54	2508118.57
30.54	2429669.31
31.54	2359427.44
32.54	2286898.45
33.54	2223763.52
34.54	2163313.5
35.54	2105424.98
36.54	2043286.52
37.54	1967678.18
38.54	1900664.32
39.54	1839481.65
40.54	1771416.04
41.54	1708788.54
42.54	1651208.45
43.54	1592274.14
44.54	1532486.8
45.54	1469841.76
46.54	1407411.7
47.54	1340524.91
48.54	1264399.52
49.54	1169862.13
50.54	1054817.24
51.54	891087.99
52.54	781758.22
53.54	621455.23
54.54	485032.18
55.54	349191.39
56.54	164635.5
57.54	62376.85
58.54	11244.41
59.54	171
)
sau <- read.csv("hipso_sau.csv")
sau <- read.csv("sau_hipso.csv")
setwd("~/ISIMIP_Lake_Sector/Hypsographics")
sau <- read.csv("sau_hipso.csv")
View(sau)
sau <- read.csv("sau_hipso.csv",header=F)
View(sau)
points(sau[,2],sau[,1], col="green",ylim = rev(range(sau[,1])))
volume <- 168555237 #m3
surface_area <- 5999766.8 #m2
max_depth <- 59.54 #m
mean_depth <- volume/surface_area
if (max_depth > floor(max_depth)){
hipsographic <- matrix(NA, floor(max_depth)+2,3)
} else {
hipsographic <- matrix(NA, floor(max_depth)+1,3)
}
z_m <- mean_depth
Vd <- 3*z_m/max_depth
b <- log10(Vd)
#if per les diferents possibilitats
Hd <- 10^(-34*b^5-18*b^4-6.3*b^3-1.9*b^2-4*b)
if (floor(max_depth)>=1){
#area across slices of the truncated cone
for (i in 0:floor(max_depth)){
hipsographic[i+1,1]=i
hipsographic[i+1,2]=surface_area*(  ( Hd^(-i/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+1,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-i/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*i/max_depth)))/(2*log(Hd)) + (Hd^-2)*i )
}
if (max_depth>floor(max_depth)){
hipsographic[i+2,1]= max_depth
hipsographic[i+2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[i+2,3]= (surface_area/((1-(Hd^-1))^2)) * (  (2*max_depth*(Hd^(-max_depth/max_depth-1)))/log(Hd) - (max_depth*(Hd^(-2*max_depth/max_depth)))/(2*log(Hd)) + (Hd^-2)* max_depth)
}
} else{
if (max_depth>floor(max_depth)){
hipsographic[2,1]= max_depth
hipsographic[2,2]= surface_area*(  ( Hd^(-max_depth/max_depth) - Hd^(-max_depth/max_depth) ) / (Hd^(-0/max_depth) - Hd^(-max_depth/max_depth) ) )^2
hipsographic[2,3]= 0
}
}
hipso_TC <-  hipsographicR_TruncatedCone(volume,surface_area,max_depth)
plot(hipsographic[,2],hipsographic[,1], col="red",ylim = rev(range(hipsographic[,1])))
lines(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,1])))
points(sau[,2],sau[,1], col="green",ylim = rev(range(sau[,1])))
plot(hipsographic[,2],hipsographic[,1], col="red",ylim = rev(range(hipsographic[,1])),log="y")
plot(hipsographic[,2],hipsographic[,1], col="red",ylim = rev(range(hipsographic[,1])),log="x")
lines(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,1])))
plot(hipsographic[,2],hipsographic[,1], type="l",col="red",ylim = rev(range(hipsographic[,1])),log="x")
lines(hipso_TC[,2],hipso_TC[,1],type="l",ylim = rev(range(hipso_TC[,1])))
points(sau[,2],sau[,1], col="green",ylim = rev(range(sau[,1])))
plot(sau[,2],hipso_TC[,2])
abline(1,1)
points(sau[,2],hipsographic[,2],col="red")
plot(sau[,2],sau[,2])
plot(sau[,2],sau[,2],type="l")
points(sau[,2],hipso_TC[,2],col="green")
points(sau[,2],hipsographic[,2],col="red")
lm(sau[,2]~hipso_TC[,2])
summary(lm(sau[,2]~hipso_TC[,2]))
summary(lm(sau[,2]~hipsographic[,2]))
serie=sau[,1]
